This calculates the length and width of the equivalent ellipse of cutouts defined appropriately within shape_definition_x_coords.csv and shape_definition_y_coords.csv. The measured cutout factors defined within measured_cutout_factors.csv are also used to calculate the equivalent ellipse.
The widths and lengths are then output into measured_cutout_factors.csv for the input into Electron Spline Modelling.ipynb.
In [ ]:
In [ ]:
In [2]:
try:
import shapely.geometry as geo
import shapely.affinity as aff
except:
!pip install shapely
import shapely.geometry as geo
import shapely.affinity as aff
try:
from descartes.patch import PolygonPatch
except:
!pip install descartes
from descartes.patch import PolygonPatch
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import pylab
%matplotlib inline
from scipy.optimize import basinhopping
from scipy.interpolate import UnivariateSpline
import os
In [2]:
In [4]:
shape_definition_x_coords = pd.read_csv("shape_definition_x_coords.csv")
shape_definition_y_coords = pd.read_csv("shape_definition_y_coords.csv")
measured_factors = pd.read_csv("measured_factors.csv")
In [ ]:
In [58]:
class equivalent_ellipse(object):
"""An equivalent ellipse given the x and y coordinates of a cutout.
"""
def __init__(self, n=5, debug=False, **kwargs):
self.debug = debug
self.cutoutXCoords = kwargs['x']
self.cutoutYCoords = kwargs['y']
self.cutout = self.shapely_cutout()
self.basinRequiredSuccess = n
self.ellipseRaw = self.ellipse_basinhopping()
self.width = self.ellipseRaw[2]
self.length = self.ellipseRaw[3]
self.ellipse = self.shapely_ellipse(self.ellipseRaw)
def shapely_cutout(self):
"""Returns the shapely cutout defined by the x and y coordinates.
"""
return geo.Polygon(np.transpose((self.cutoutXCoords,self.cutoutYCoords)))
def shapely_ellipse(self,ellipseRaw):
"""Given raw ellipse values create a shapely ellipse.
"""
xPosition = ellipseRaw[0]
yPosition = ellipseRaw[1]
width = ellipseRaw[2]
length = ellipseRaw[3]
rotation = ellipseRaw[4]
unitCircle = geo.Point(0,0).buffer(1)
ellipse = aff.scale(unitCircle, xfact=width/2, yfact=length/2) # Stretched
ellipse = aff.translate(ellipse, xoff=xPosition, yoff=yPosition) # Translated
ellipse = aff.rotate(ellipse, rotation) # Rotated
return ellipse
def minimise_function(self, ellipseRaw):
"""Returns the sum of area differences between the an ellipse and the given cutout.
"""
ellipse = self.shapely_ellipse(ellipseRaw)
return ellipse.difference(self.cutout).area + self.cutout.difference(ellipse).area
def ellipse_basinhopping(self):
"""Fitting the ellipse to the cutout via scipy.optimize.basinhopping.
"""
self.functionReturns = np.empty(self.basinRequiredSuccess)
self.functionReturns[:] = np.nan
self.numSuccess = 0
minimizerConfig = {"method": 'BFGS'}
initial_input = np.array([0,0,3,4,0])
basinhoppingOutput = basinhopping(self.minimise_function,
initial_input,
niter=1000,
minimizer_kwargs=minimizerConfig,
take_step=self.step_function,
callback=self.callback_function)
return basinhoppingOutput.x
def step_function(self,optimiserInput):
"""Step function used by self.ellipse_basinhopping.
"""
optimiserInput[0] += np.random.normal(scale=1.5) # x-position
optimiserInput[1] += np.random.normal(scale=1.5) # y-position
optimiserInput[2] += np.random.normal(scale=3) # width
optimiserInput[3] += np.random.normal(scale=4) # length
optimiserInput[4] += np.random.normal(scale=90) # rotation
return optimiserInput
def callback_function(self, optimiserOutput, minimiseFunctionOutput, minimiseAccepted):
"""Callback function used by self.ellipse_basinhopping.
"""
if self.debug:
print(optimiserOutput)
print(minimiseFunctionOutput)
print(minimiseAccepted)
print(" ")
if minimiseAccepted:
if self.numSuccess == 0:
# First result
self.functionReturns[0] = minimiseFunctionOutput
self.numSuccess = 1
elif minimiseFunctionOutput >= np.nanmin(self.functionReturns) + 0.0001:
# Reject result
0
elif minimiseFunctionOutput >= np.nanmin(self.functionReturns) - 0.0001:
# Agreeing result
self.functionReturns[self.numSuccess] = minimiseFunctionOutput
self.numSuccess += 1
elif minimiseFunctionOutput < np.nanmin(self.functionReturns) - 0.0001:
# New result
self.functionReturns[0] = minimiseFunctionOutput
self.numSuccess = 1
if self.numSuccess >= self.basinRequiredSuccess:
return True
In [ ]:
In [4]:
In [ ]:
In [ ]:
In [ ]:
In [48]:
test_x = shape_definition_x_coords["Simon #3 cutout"].values
test_x = test_x[~np.isnan(test_x)]
test_x
Out[48]:
In [49]:
test_y = shape_definition_y_coords["Simon #3 cutout"].values
test_y = test_y[~np.isnan(test_y)]
test_y
Out[49]:
In [59]:
test = equivalent_ellipse(n=2,debug=True,x=test_x,y=test_y)
In [66]:
plt.plot(test.cutout.exterior.xy[0],test.cutout.exterior.xy[1])
plt.plot(test.ellipse.exterior.xy[0],test.ellipse.exterior.xy[1])
plt.axis("equal")
Out[66]:
In [61]:
test.ellipse
Out[61]:
In [62]:
test.cutout
Out[62]:
In [ ]:
In [12]:
a = np.random.normal(size=10)
b = np.random.normal(size=10)
In [13]:
a
Out[13]:
In [14]:
b
Out[14]:
In [31]:
list(tuple(np.transpose([a,b])))
Out[31]:
In [39]:
geo.Polygon(np.transpose((a,b)))
Out[39]:
In [ ]: